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Abstract 

We present a detailed study of rectilinear shear deformation in 
the framework of orthotropic nonlinear elasticity, under Dirichlet and 
mixed-boundary conditions. We take a slab made of a soft matrix, 
reinforced with two families of extensible fibers. We consider the case 
where the shear occurs along the bissectrix of the angle between the 
two privileged directions aligned with the fibers. We show that if the 
two families of parallel fibers are mechanically equivalent, then only 
smooth solutions are possible, whereas if the mechanical differences 
among the two families of fibers is pronounced, then strain singu- 
larities may develop. We determine the precise conditions for the 
existence of such singular solutions for the standard reinforcing or- 
thotropic model. We then extend our findings to some orthotropic 
models of interest in biomechanical applications, and we discuss the 
possible relevance of the singular solutions to biomechanics. 



1 Introduction 



Biological soft tissues exhibit complex mechanical behaviors, which are not 
easily accounted for by classic elastomeric constitutive models. The extension 
of the mathematical models of nonlinear elasticity from rubber to soft tissues 
continues to be a challenging area in theoretical biomechanics. For example, 
the presence of oriented collagen fiber bundles in blood vessels calls for the 
consideration of anisotropy in the mathematical modeling of the mechanics 
of arterial tissues, but the mathematical theory of nonlinear hyperelastic 
anisotropic materials is not as developed as the theory of isotropic nonlinear 
elasticity. A consultation of the eminent book by Antman [1] shows that, 
in recent years, there has been very few additions to the classical works 
of Adkins |2] or of Ericksen and Rivlin [3] with respect to the solution of 
boundary-value problems in nonlinear anisotropic elasticity. 

It is well known that certain radial anisotropies in linear and non-linear 
elasticity problems can give rise to stress singularities which are absent in the 
corresponding isotropic version of these problems. Lekhnitskii [J| was per- 
haps the first to observe this peculiarity, by studying a circular orthotropic 
plate compressed by a uniformly distributed force along its external bound- 
ary. Antman and coworkers (see for example |5]) extended in some sense this 
analysis to radially symmetric equilibrium states of anisotropic nonlinearly 
elastic bodies. Another example of this extension to nonlinear elasticity is 
found in the paper by Kassianadis et al. [U] on the finite azimuthal shear of 
transversely isotropic materials. 

Merodio et al. [7j investigated a simple model for a nonlinear, trans- 
versely isotropic, elastic solid and discovered a new kind of singular behavior, 
not present in isotropic materials. It occurs for the inhomogeneous rectilin- 
ear shear of an incompressible elastic slab reinforced by a family of parallel 
fibers. They show that, depending on the reinforcement strength and on the 
fiber orientation with respect to the shearing direction, weak solutions for 
this simple boundary value problem may be expected. These solutions are 
associated with fiber kinking and loss of ellipticity of the field equations. The 
deformation field is continuous, but it suffers a jump in the first derivative 
and a blow-up for the second derivative. Therefore the stress field suffers a 
discontinuity of first kind, a phenomenon clearly associated with mechani- 
cal instabilities. It also puts into question the applicability of finite element 
methods to nonlinear anisotropic elasticity, because the obtention of numer- 
ical solutions to the governing equations requires the calculation of second-, 
fourth-, and sometimes higher-order derivatives. In biomechanical applica- 
tions of the constitutive models of arterial walls, the appearance of stress 
singularities is an important mathematical aspect of the theory, because it 
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may be associated with some pathological states of the tissues (such as the 
bursting of an aneurysm). 

The aim of the present paper is to extend the results of Merodio et al. 
[7] from transverse isotropy (one family of parallel fibers) to orthotropy (two 
families of parallel fibers), often encountered in biological soft tissues. We 
note that Fosdick and Royer-Carfagni |8j show that Lekhnitskii's classical so- 
lution predicts the interpenetration of material regions, an unacceptable de- 
formation behavior in the classical theory of elasticity. However the solutions 
proposed in [3| and here are isochoric and thus satisfy the local injectivity 
requirement. 

The paper is organized as follows. In the next section we write down 
the governing equations and boundary conditions, and we discuss the basic 
mathematical issues at play. Section 3 is devoted to one of the simplest model 
of nonlinear orthotropic elastic materials (the standard reinforcing model), 
obtained by adding the classical neo-Hookean strain energy density to two 
terms that take into account the reinforcements along the fiber directions. 
These latter terms are quadratic in the squared extension along the fibers. 
We solve the problem of inhomogeneous rectilinear shear along the bissec- 
trix to the fibers, first for Dirichlet boundary conditions and next for mixed 
boundary conditions. We also provide an energy analysis of the solutions. 
In Section 4 we consider a more advanced constitutive model of the biome- 
chanics literature, proposed by Holzapfel et al [9], where the reinforcement 
terms in the strain-energy density are exponential, in order to account for a 
strong stiffening effect (the artery model). 

The results suggest that orthotropic fiber reinforcement is quite efficient 
at cancelling the singularities and the shear discontinuities encountered in 
transversally isotropic fiber reinforcement. Indeed we recover the main fea- 
tures discovered by Merodio et al. [7j (jump in the shear, blow-up of the 
second derivative of the displacement) but under the condition that one 
family of fibers is much stiffer than the other. For the standard reinforc- 
ing model, one stiffness modulus must be at least 9.9 times larger than the 
other; for the artery model, the stiffnesses ratio is even higher, due to ex- 
ponential terms. In general, the families of parallel collagen fibers found in 
arteries are determined experimentally to be mechanically equivalent, sug- 
gesting that singularities do not develop, at least in physiological conditions, 
for the rectilinear shear of arteries. 
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2 Basic equations 



We consider a composite incompressible slab with thickness L, made of an 
isotropic matrix reinforced with two families of parallel extensible fibers (the 
fibers are all orthogonal to the boundaries of the solid.) In the undeformed 
configuration, we call {Xi, X2, X^) a set of Cartesian coordinates such that 
the solid is located in the ^ X3 ^ L region. We denote by Ei, E2, 
the orthogonal unit vectors defining the Lagrangian (reference) axes, aligned 
with the Xi, X2, X3 directions, respectively. 

When the solid is sheared in the direction of Ei, the particle initially at 
X moves to its current position x. We call F = dx/dX the associated defor- 
mation gradient tensor, and B = F^F the left Cauchy-Green strain tensor. 
We then call (xi, X2, x^) the Cartesian coordinates, aligned with (Xi, X2, X3), 
corresponding to the current position x. In the current configuration, the 
basis vectors are ei, 62, 63, and here they are such that = Ei {i = 1, 2, 3). 
The deformation is given in all generality by 

Xi =Xi+L/(X3/L), X2=X2, X3 = X3, (2.1) 

where / is a yet unknown function oi t] = X3/L only. The amount of shear 



is /' = df /dr]. The deformation (2.1 ) is a simple shear when /' is a constant; 
otherwise it is a rectilinear inhomogeneous shear. The direction of shear is 
that of Ci = El and the plane of shear is that of (ei = i^i, 62 = -E2). 
We find in turn that 

F = I + fei ®E3, B = I + /(ei ® 63 + Ci ® 63) + (ffei ® Ci. (2.2) 

The first principal isotropic strain invariant Ii = ti B is given here by 

/l = 3+(r)^ (2.3) 

and the second principal isotropic strain invariant, I2 = [If — tr (B'^)]/2, is 
also equal to 3 + {f'Y- 

We call $ (\l/, respectively) the angle between the direction of one family 
of parallel fibers (the other family, respectively) and the direction of shear 
Xi. In other words, the unit vectors M and N (say) in the two preferred 
fiber directions have components 

M = cos^Ei + sin^Ea, N = cos'^Ei + sin'^Es, (2.4) 

and they are transformed into m = FM and n = Fn in the current con- 
figuration, 

m = (cos $ + /' sin $)ei + sin $63, n = (cos ^ + f sin \I')ei + sin $63. 

(2.5) 
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Figure 1 : Two unit squares lying in the transverse section of a slab reinforced with 
two families of fibers (thin lines) and subject to a simple shear of amount 0.5 along 
the bissectrix of the angle between the two families. In the reference configuration, 
one family of fibers is aligned with the unit vector M making an angle <I> = 60° 
with the Xi-axis; the other family is aligned with N, at an angle 120° with the 
Xi-axis. In the current configuration, they are along m and n, respectively. 

In the remainder of the paper, we restrict our attention to the special case 
where the material is sheared along a bisectrix of the angle between the two 
families. Generahty is lost with this approach, but it has the merit of keeping 
low the number of geometric parameters; we also argue that it still captures 
some salient features of sheared soft tissues with two preferred directions. 

Hence from now on, \1/ = vr — $ and the angle between the two preferred 
directions is tt — 2$. In other words, the unit vectors M and N in the 
preferred fiber directions have components 

M = cos$^;i + sin<l>^;3. A/" = -cos<l>^;i + sin^E^s, (2.6) 

in the reference configuration, and they are transformed into 

m = (cos $ + /' sin $)ei + sin $63, n = (— cos $ + /' sin $)ei + sin $63, 

(2.7) 

in the current configuration; Figure [T] is a visualization of the situation in the 
case of a simple (homogeneous) shear of amount 0.5 and an angle $ = 60°. 
Note that because the reinforcements are not directional, we may without 
loss of generality restrict ourselves to the range < $ < vr. 

We now introduce the anisotropic invariants 1^ = m • m and = 
Fm-Fm; in particular we find 

/4 = l + /'sin2<l> + (/')^sin2<l>. (2.8) 

Recall that I4 is the squared stretch in the fiber direction (Spencer [10]). In 
particular, if I4 ^ 1 then the fibers aligned with m are in extension, and 
if /4 ^ 1 then they are in compression. Clearly here, when ^ $ ^ 7r/2, 
the quantity /4 — 1 is always positive and the fibers aligned with m are in 
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extension. On the other hand, when 7r/2 < $ < tt, there always exist a 
certain amount of shear (exphcitly, —2/ tan$) above which these fibers are 
in compression. 

The other anisotropic invariants are Iq = n • n, = Fn-Fn, and Ig = 
m • n. Here we find that 

Je = 1 - /'sin2$ + (/')^sin2$. (2.9) 

In general, the strain-energy density W of a hyperelastic incompressible 
solid reinforced with one or two families of parallel extensible fibers depends 
on two isotropic deformation invariants: Ji and I2, and on the five anisotropic 
deformation invariants (TUl E]: I4, . . . , Is- Henceforth we make the assump- 
tion that W is the sum of an isotropic part and an anisotropic part. For the 
isotropic part, modelling the properties of the elastin matrix, we take the 
neo-Hookean strain-energy density, with constant shear modulus fi. For the 
anisotropic part, modelling the properties of the extensible collagen fibers, 
we take the sum of a function of J4 only and a function of Jg only, say 
F(/4) + G{Iq). Hence we restrict our attention to those solids with strain 
energy density 

W = fi{h-3)/2 + Fih) + Gih). (2.10) 

Now the Cauchy stress tensor cr derived from this strain energy function 
is (see e.g. Ogden, 1984), 

cr = -pi + 12B + 2F'{h)m ® m + 2G'{I(i)n ® n, (2.11) 

where p is a Lagrange multiplier introduced by the constraint of incompress- 
ibility, and F' = dF/d/4, G' = dG/d/4. 

Because shear is a plane strain deformation, and because the fibers lie in 
the plane of shear, it is a simple matter to find the directions of principal 
stresses. One is normal to the plane of shear, and the two others are in the 
(ci, 63) plane, at the angles (p and ip + TT/2 from the direction of shear, where 
ip g]0, 7r/4] is defined by 

tan 2(f = 2(ei • cre^,) / (ei ■ crci — 63 ■ cre^). (2-12) 

Here we find that 

ei ■ (763 = /if + 2F{h)m,m3 + 2G'{h)n^n^, (2.13) 



where mi = m • and rii = n • are found from (2.7). 
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The equilibrium equations, div cr = (in the abscence of body forces) 
reduce to 



1^ = [fif + 2F'ih)mm3 + 2G'ih)n,ns] , (2.14) 
ox I ax 3 

0, 



dx2 

|^ = ^[/i/' + 2F'(/4K + 2GUK], 

where the expressions in brackets are independent of xi and X2- It follows 
that 

p = p(xi, X2) = Cosi + 2F'{h)ml + 2G'{h)nl + D, (2.15) 

where Cq, D are arbitrary constants of integration, is a suitable pressure field. 
A single governing equation remains to be solved for the shear deformation, 
namely 

^ [/i/ + 2F'(/4)mim3 + 2G\h)n^n3] = C^L. (2.16) 

We consider two specific boundary value problems (BVPs). In the refer- 
ence configuration, the slab of thickness L in the X3 direction and of infinite 
dimensions in the other directions is bonded to two infinite rigid plates lo- 
cated at X3 = and X3 = L. A constant pressure gradient is applied in 
the Xi direction and drives the deformation of the slab. The overall goal 



of our investigation is to solve (2.16) subject (i) to the Dirichlet boundary 



conditions: /(O) = 0, /(I) = 0, and (ii) to the mixed boundary conditions: 
/(O) = 0, /'(I) = Ki, where Ki is a prescribed constant. In Case (i), we have 
a classical two-point boundary value problem which, for an isotropic medium, 
may be reduced to a Cauchy problem by using symmetry considerations. In 
our anisotropic case it may happen that, once the second-order differential 



equation (2.16) is rewritten in normal form, the corresponding right hand- 
side is neither continuous nor Lipschitzian with respect to /'. Then standard 
methods for the study of the existence and uniqueness of the solution may not 
apply any longer; moreover, the solution may develop singularities. In Case 
(ii), the BVP is simpler to solve, but it is also possible to have non-smooth 
solutions. We point out that enforcing the boundary condition /'(I) = Ki is 
equivalent to prescribing the shear stress T12 on the upper face of the slab. 

For transversally isotropic materials. Case (i) has been studied by Mero- 
dio et al [7\ and a mixed-BVP similar to Case (ii) has been considered for 
azimuthal shear by Kassianadis et al [6]. 
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3 The standard reinforcing model 
3.1 Normal form of the BVP 

The standard reinforcing model for solids with two family of fibers is a special 
case of (2.10). Its strain energy density is 

W = ^(/i - 3)/2 + fiE^ih - 1)V4 + /i^2(/6 - 1)V4, (3.1) 

where fiEi and fiE2 are the extensional moduli in the fiber directions. The 
BVPs based on (2.16) are now 

^ [/' + E,{h - l)mm3 + E2{h - Ihms] = CoL/^i; (3.2) 

with the boundary conditions (i): /(O) = 0, /(I) = 1 and (ii): /(O) = 
0, /'(I) = Ki. We begin our study with the Dirichlet boundary conditions, 
Case (i). 

The differential equation may be rewritten as 

^ {/' + 7/' sin^ *[2 cos^ $ + 3/3(/') sin $ cos $ + {ff sin^ 0] } = CoL//i, 

(3.3) 

where we introduced the dimensionless material constants 7 and /5, defined 
as 

7 = ^1 + ^2, (3 = {E,- E2)/iE^ + E2). (3.4) 

The quantity 7 gives a measure of the collagen/elastin strength ratio, and 
the quantity /3 gives a measure of the orthotropy. If 7 = 0, then the material 
is isotropic. If /3 = ±1, then either ii^i = or £^2 = and the solid is 
transversally isotropic (there is only one active family of parallel fibers); if 
/3 = 0, then Ei = E2 and the two families of fibers are said to be mechanically 
equivalent. 

In its normal form, the BVP Case (i) reads 

where a = CqL/h is a dimensionless measure of the pressure gradient and 
where the denominator D is defined as 

D{f\^) = l + 7sin2$[2cos2$ + 6/3cos$sin$(/) + 3sin2$(/')2]. (3.6) 

First we note that when /3 = 1, the whole analysis is consistent with 
that of Merodio et al. [Tj for a transversally isotropic slab. Also, when 
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/3 = 0, the governing equation coincides with that obtained for the rectihnear 
inhomogeneous shear of the isotropic sohd slab with strain energy density 

= (/i + 27 sin^ $ cos^ $)(/!- 3)/2 + (7sin^ - 3f/A. It follows that 
when the two families of fibers are mechanically equivalent, only smooth 
solutions exist and no singularity may develop. 

Next we take /3 7^ and notice that D is a quadratic in the amount of 
shear /'. If its discriminant is negative, then no singularity may develop. 
The denominator D has real roots when 

(3^2 -2)7sin2 2$ -4 > 0. (3.7) 

Therefore a necessary condition for the appearance of singularities is that 

> 2/3. (3.8) 

Assume that the fibers along M are stiffer than those along N. Then Ei > 
E2 and this inequality means that E1/E2 > 5 + 2^ ~ 9.9. Hence we are 
certain that singularities do not develop when the fibers along M are less 
than 9.9 times stiffer than the fibers along N . 



3.2 Orthogonal fibers: $ = 7r/4 

Here we focus on the special case where one family of fibers is orthogonal to 



the other family ($ = 7r/4). Then the denominator D in (3.6) reduces to 



D{f, 7r/4) = 1 + (7/4) [2 + 6/3(/') + 3(/)^] . (3.9) 

Clearly, whether /" develops singularities or not depends among other things 
on the sign of the quantity (3/3^ — 2)7 — 4. Figure |2] displays on the left the 
curve where this quantity is zero in the (/3,7) plane. When it is negative, 
the existence and uniqueness of a smooth solution are guaranteed by general 
theorems and standard numerical procedures of integration can be imple- 
mented. For instance, we take (3 = 0.5, 7 = 3.0, and a = 1.0,5.0,10.0 in 
turn, and obtain the displacements displayed on the right of Figure [2| using 
the finite difference method implemented into Maple. 

When (3/3^-2)7-4 ^ 0, there is a chance that singularities may develop 
within the thickness of the slab and we now investigate this possibility. First 
we consider the case where this discriminant is equal to zero, when 7 = 
4/(3/32-2). Then 

D = ^P^(f' + ^)"^ (3.10) 



and integrating (3.5) once gives 

(/' + /3)3 = a{3p' - 2)(r/ - r^o) + /3^ (3.11) 
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Figure 2: On the left: curve in the (/S, 7) plane separating the region where only 
smooth exist ((3/3^ — 2)7 — 4 < 0) from the region where singularities might develop 
for the second derivative of the displacement, in the case where the two families of 
fibers are orthogonal. On the right: an example of a completely smooth solution, 
obtained for /3 = 0.5, 7 = 3.0 and for several values of the pressure gradient, as 
measured by a. 



where ?7o G (0, 1) is a point in the thickness of the slab where /' = (its 
existence is ensured by the continuity and differentiability of /, coupled to 
the boundary conditions /(O) = /(I) = 0). Solving for /' gives 



(3.12) 



where the sign depends on the sign of the radical. Integrating further, and 
imposing /(O) = 0, we obtain 



fiv) 



4a(3/32 - 2) 



-[P^-a{3P'-2)r]of'}-Pr]. (3.13) 



To solve the BVP entirely, it remains to determine r]o G (0, 1). It is fixed by 
the second boundary condition: /(I) = 0, i.e. it is a solution to the equation 



[^^ + «(3/32 - 2)(1 - r^o)]'^' - [P' - «(3/3' - 2)r]of' = Aa/3{3P^ - 2)/3. 

(3.14) 



4/3 



Now, collecting (3.5), (3.10), (3.11), we see that /" blows up at = rjs given 
by 



a{3P'-2){7]s-Vo)+P' 



0. 



(3.15) 
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Figure 3: Curves in the {r]Q, a) plane giving the locus of the constant of integration 
7]Q corresponding to a given level of pressure gradient, as measured by a, in the case 
where the two families of fibers are orthogonal and some fibers are stiff enough 
to guarantee the appearance of singularities. The curves are shown in the first 
quarter of the plane, for /3 = 0.875,0.9,0.925,0.95,1.0. They are antisymmetric 
with respect to the point (0.5,0) (second part not shown here). The curves are 
limited to the right by a vertical line at rjQ = 27/64 = 0.421875. 



The final condition to impose for this singularity is that it occurs within the 
thickness of the slab: ^ 775 ^ 1, i.e. 



a(3/32 



^ 1. 



(3.16) 



On Figure |3] we graph the curves defining the pairs {rjoja) such that the 
second boundary condition (3.14 ) is satisfied, for several values of (3. We limit 



the display to the range r/o < 0.5 because the curves are antisymmetric with 
respect to the point (0.5,0). For visual reasons, the upper bound is taken 
as ajnax = 40.0. The other limit of each curve is imposed by the inequal- 
ities (3.16), specifically here the lower one. The corresponding transitional 



behavior is dictated by the equality 

VO — (^70) trans 



a(3/32-2)' 



(3.17) 



When this holds, the corresponding transitional level of pressure gradient is 



found from (3.14) as 



a 



(Cl) trans — 



64/3^ 
27(3/32 - 



(3.18) 
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Figure 4: Plots of the displacement and of its first and second derivatives through 
the slab thickness in the case where the two families of fibers are orthogonal and 
singularities develop. Here (3 = 0.875 (large difference in fiber stiffness) and a 
(giving a measure of the pressure gradient) and ryo (constant of integration) are 
chosen so that the second derivative is discontinuous. Thin curves: a = 5.3489, 
rjQ = 0.421875); Medium thickness curves: a = 10.0, ryo = 0.33120; Thick curves: 
a = 15.0, r/o = 0.30871. 



Substituting back above gives 

(//o)trans = 27/64 = 0.421875. (3.19) 

Hence all the curves stop at the vertical barrier rj^ = 0.421875, irrespective 
of the value of (3. 

For all values of a and t]o such that the point {a, r]o) belongs to one of these 
curves, a singularity develops within the thickness for the second derivative 
of the displacement. For instance at the points ( (a)trans, ('7o)trans), the exact 



solution (3.13) and its derivatives reduce to 



/(r^) = (3 {v'/' - v) , f'iv) = /3 [(4/3)r^^/V3 - l] , fiv) = (4/9)/3r^/^ 

(3.20) 

and /" clearly blows up on the r] = face of the slab. For Figure |4] we take 
P = 0.875 and three points on the corresponding curve of Figure |3] , namely 
(a = 5.3489, r/o = (r?o)trans), (a = 10.0, r^o = 0.33120), and (a = 15.0, r/o = 
0.30871). For the first combination, /" blows up on the slab face 7] = 0; for 
the second and third combinations, it blows up within the thickness of the 
slab. 

Next we consider the case where the discriminant of the /' quadratic in 



(3.9) is positive: (3/3^ — 2)7 — 4 > 0, and focus now on finding singularities 
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9(x) 

Figure 5: (a) Plot of the cubic establishing a relationship between the amount of 
shear and the thickness. When gi ^ ij ^ §2, there are three possible values of /' 
for each r]. The vertical dashed lines are at rj = 51,52; the vertical full line defines 
two regions of equal area between gi and 52- (b) Maxell rule convention: jump 
in the amount of shear at the thickness giving equal areas, (c) Maximum delay 
convention: jump in the amount of shear at rj = gi. Here the Dirichlet BVP is 
solved for j3 = 0.95, 7 = 10.0, a = 10.0. 



for /', the amount of shear. Integrate the BVP (3.5), (3.9) once to get 



where ?7o is a constant of integration, and g is the following cubic, 

/ N f, , ^\ , 3/37 2 , 7 3 

g{x) := \ l + -jx+ —X + -x^ 
with a local maximum (resp. minimum) at Xi (resp. X2) defined as 



(3.21) 



(3.22) 



3^1,2 



-It 



'(3/32-2)7-4 



For gi = g{xi) and g2 = g{x2), we find 
1 



37 



7 



(3.23) 



g,,2 = -g [2(2 + /3) + 7(1 - /3)(3/3 - 2)] ± ^ 



(3/3^-2)7-4 ' 
37 



3/2 



. (3.24) 



Clearly, a systematic procedure to establish a one-to-one correspondence 
between x and g everywhere (or equivalently, between /' and 77) runs into 
difficulties in the interval xi ^ x ^ X3, where 0:3 is the root of g{x) = gi 
other than xi, see Figure [5|^a). To address this problem, we take the stance 
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that /' jumps from a low value to a higher one. In order to jump follow- 
ing the absolute minima of the energy, we consider in turn the Maxwell rule 
convention of equal area, see Figure |5](b), and the Maximum delay conven- 
tion, see Figure |5](c). We propose to track these two possible solutions by 
a suitable numerical approach. This hands-on approach is required because 
usual numerical methods sometimes fail in finding good approximations. In 
fact, commercial code solvers issue a warning here about possible failure in 
the numerical convergence and are unable to provide a satisfactory solution 
in this region. Note that the non-monotonous behavior does not necessar- 
ily occur within the slab thickness and that some parameter values allow a 
monotonous variation of /', devoid of jumps (such is for instance the case 
when Qi ^ 1). We focus on those parameter values which do give a jump 
inside the slab. 

The main difficulty in solving the Dirichlet BVP is that we do not have 
an analytical access to the value of the integration constant r]o in (3.21). We 
tackle the Dirichlet BVP by a shooting method, combined with the bisection 
method, in the following manner. 

We take r^Q^-* (say) as an initial guess for rjQ. Then let Kq^^ = /'(O); it is the 
real root to the cubic g{K^'') = — ar/g^''. It is now possible to reformulate the 
BVP as an Initial Value Problem (IVP), which we solve numerically on two 
subintervals of [0, 1]. That process is detailled later, in the simpler case of the 
mixed BVP. It gives ri^^\ the thickness where the jump takes place, and also 
f{ri) numerically. Finally we compute /(I) and measure how different it is 
from the second boundary condition /(I) = 0: if |/(1)| ^ tol is not satisfied, 
for a prescribed numerical tolerance "tol", then we adjust the approximate 
value of rjf) from r/g^^ to rjl^^ and so on, from r]^'^"* to ?7o until the criterion 
of convergence is reached, following the indications given by the bissection 
method. In the process we also get access to r^^ , a numerical approximation 
of the singularity point 775. 

In Figure [6} we report the numerical solutions for 7 = 10.0, a = 10.0, 
and in turn, (3 = 1.0 and /3 = 0.95, with tol=le-6 as the tolerance for 
the stopping criterium in the bissection method. The values identified by 
the bisection method for the integration constants are as follows. When 
P = 1.0 (transverse isotropy), we find rjo = 0.2423 for the Maxwell rule 
solution and 7]o = 0.21725 for the Maximum delay solution; when P = 0.95 
(orthotropy) , we find r/o = 0.13818 for the Maxwell rule solution and tiq = 
0.05014 for the Maximum delay solution. It is worth noting that the two 
kinds of solutions not only jump at different singular points, but also present 
different slopes, before and after the singular points. Moreover, we checked 
that the solutions obtained for /3 = 1.0 (transverse isotropy) are consistent 
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Figure 6: Maxwell rule solutions and Maximum delay solutions for /3 = 1.0 (trans- 
verse isotropy) and for f3 = 0.95 (orthotropy) , with zero Dirichlet boundary con- 
ditions. Here 7 = 10.0, a = 10.0. 



with those obtained by Merodio et al. [7], using a different numerical method, 
based on a quadrature approach. 

We now consider the mixed BVP, Case (ii), 

which turns out to be simpler to analyze and to solve numerically than the 
Dirichlet BVP. 

The main features uncovered in the previous analysis still apply. Hence 
the uniqueness of the solution is not guaranteed for all parameter values, 
because the energy can have two minima and is in general not a convex func- 
tion, leading to a jump in the derivative of the displacement. The analysis 
for the mixed boundary conditions is almost identical to that of the Dirichlet 
boundary conditions, with the difference that it is now possible to identify a 
priori the location of the singularity. 

In order to jump following the absolute minima of the energy, again we 
consider in turn the Maxwell rule convention of equal area, and the Maximum 
delay convention, because commercial solvers also fail here. We track these 
solutions by transforming the mixed BVP into a second-order initial value 
problem (IVP), as follows. 



Starting from the first integral (3.21 )-(3.22) of equation (3.3), we find 



from the second boundary condition /'(I) = Ki that 

r/o = l-(7(i^i)/a. (3.26) 
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Then let Kq = /'(O); it is the real root to the cubic 



g{Ko) = -ar]o = g{Ki) - a. (3.27) 

Now we can reformulate the BVP as an IVP, which we solve numerically in 
two steps. First on the subinterval [O,//^], with initial conditions: /(O) = 0, 
/'(O) = Kq] we call fs and Ks the computed values of / and f at r] = rjs, 
the slab thickness where the jump takes place. Next we solve numerically 
the second part of the IVP, this time on the subinterval [r]s, 1], with initial 
values: /(r/s) = fs, f'ijis) = Ks- To compute the value of rjs, the singularity 
thickness, we proceed as follows. 

In the case of the Maxwell rule convention of equal area, the singularity 
occurs at the inflection point of the function g. Solving g"{Ks) = gives 
Ks = —(3 and then, g{Ks) = /3[7(/3^ ~ l)/2 — !]• Then r]s is found by solving 
the equation 

giKs)=giK,) + air]s-l). (3.28) 
In the case of the Maximum delay convention, the singularity occurs at 



the local maximum of the function g. Hence Ks = Xi given by (3.23); then 



g{Ks) = gi given by (3.24), and rfs is found from (3.28) for 775. 

Figure [7] shows the numerical solutions obtained with this numerical tech- 
nique for the values 7 = 10.0, a = 10.0, Ki = 0.5, and in turn, /3 = 1.0 
(transverse isotropy) and /3 = 0.95 (orthotropy) . In the figure on the left, we 
report the numerical approximation for /(r/) and in the figure on the right, 
the approximations for the amount of shear /'(?]), clearly showing that the 
jumps of the derivatives occur at different singular points. For that example, 
we find rjo = 0.48125 when f3 = 1.0 and r]o = 0.44375 when /3 = 0.95. 



3.3 Non-orthogonal fibers: $ 7^ 7r/4 

Extending the results and techniques developed at $ = 7r/4 to the case 
$ 7^ 7r/4 (non-orthogonal fibers) poses no particular problem. Rather than 
detailing the process, we refer the reader to the paper by Merodio et ai, [7j 
where the extension is done in the case /3 = 1 (transverse isotropy). 



4 Orthotropic biomechanical model 

We now investigate briefly whether the analysis conducted for the stan- 
dard reinforcing model can be extended to a strain energy density often 
encountered in the biomechanics literature, namely the model proposed by 
Holzapfel et al. [9J to describe the behavior of an orthotropic artery, and 
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Figure 7: Plots of the displacement and of its first derivative through the slab 
thickness in the case where the two families of fibers are orthogonal. Numerical 
solutions for the mixed BVP obtained by following the Maxwell rule convention 
(round dots plots) and the Maximum delay convention (square dots plots), for 
(3 = 1.0 (transverse isotropy) and /3 = 0.95 (orthotropy), and 7 = 10.0, a = 10.0, 
Ki = /'(I) = 0.5. 

widely used since, for instance to model porcine aortic tissue, passive basilar 
artery, cornea, etc. We present it in the form 



iy=^(/i-3) + ^{exp [h{h-ir]-l} 

fiE: 



+ ^{exp[fc,(/e-l)Vl}> (4-1) 



where fci, k2 are dimensionless constants. Equation (2.16) is then rewritten 
as 



drj 



{/' + E,{h - 1) exp [hih - 1)'] mm3 



+E2{Ie - l)exp [hih - If] n^n^} = CoL/fi. (4.2) 



The BVP can be put in the form (3.5), where now 



DU", $) = 1 + sin $ {EiTiif, $) exp [h{h - if] 

+E2T2if',^)exp[k2ih-iy]}, (4.3) 
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and the functions Fi and r2 are defined as 



Ti = 2kif'^ (/' sin $ + cos $) (/' sin^ $ + sin 2$) ^ (2/ sin^ $ + sin 2$) 

+ 3(/')^ sin^ $ + 3/ sin $ sin 2$ + cos $ sin 2$, (4.4) 

and 

12 = 2^2/2 (/' sin $ - cos $) (/' sin^ $ - sin 2$) ^ (2/ sin^ $ - sin 2$) 

- 3{f'y sin^ $ + 3/' sin $ sin 2$ + cos $ sin 2$. (4.5) 

To simplify the algebra we restrict our discussion to the special case where 
the families of fibers are at right angle, $ = it /A. Our objective is to find out 
if there exist special values /' = (/')* say, such that D{{f')*, 7r/4) = 0. Now 
Fi and r2 reduce to 

ri(/', rr/A) = ^ {r (/' + iy{f + 2fk, + (3^ + 6/ + 2)} , 

r2(/', vr/4) = ^ [fif - Ifif - 2fh + (3^ - 6/' + 2)} . (4.6) 



In the biomechanical applications of the model (4.1), it is often assumed 
that the two families of fibers are mechanically equivalent, so that Ei = E2 
and ki = ^2- In that case, some long but simple computations show that 
/' = is a minimum for the function D. Because D(0,7r/4) = 1 + i^i 7^ 0, 
we conclude that singularities may not develop (This result may be extended 
to any angle $ quite easily). 

When El 7^ E2, things are more complex. For instance, consider the 
values of D{f', 7r/4) when /' = —1, 0, 1 in turn: 

D{-1, 7r/4) = 1 - Eie"'/^ + (36^2 + ll)(^2/4)e='*^^/^ 
D(0,7r/4) = 1 + (El + E2)/2, 

D{1, 7r/4) = 1 + (36A;i + ll){Ei/4:)e^'''/^ - E^e^^'^ (4.7) 
Therefore, because -D(0,7r/4) > 0, it is sufficient to choose 

9,Qk + 11 

exp(A;i/4)Ei > 1 + exp(9A;2/4)E2 (4.8) 

to obtain the existence of at least one (/')* such that D{{f')*, 7r/4) = 0. This 
inequality suggests that singularities occur only for huge differences between 
the fiber stiffnesses, and are unlikely to be observed at all for realistic values 
of the parameters. Take for example the case where ki = k2 = k (say). Then 
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Figure 8: Numerical solution to the zero Dirichlet BVP for the aartery model. 
Here ki = k2 = 0.1 and /3 = 1.0 (transverse isotropy), 0.9, 0.875, 0.82 (orthotropy) . 
Other parameters: Ei + E2 = 10.0, CqL/h = 10.0. 



P defined in (3.4) gives a measure of the orthotropy: when /3 = 1, the sohd 
is reinforced with one family of parallel fibers, and when /3 < 1, there are 
two families of parallel fibers at play. To generate the graphs in Figure |8| 
we take = 0.1 and /3 = 1.0, 0.9, 0.875, and 0.82 in turn. At /3 = 1.0, the 
shear variations are pronounced but regular. As soon as /3 < 1 (two families 
of fibers), the shear variations are quickly smoothed down, highlighting the 
stabilizing effect of orthotropy. 

We now evoke some possible applications of our results to biomechanics. 
Indeed, we know that arterial tissue adapts to physiological and pathological 
stimuli though rearrangement of the microstructure. Arterial remodeling is 
induced by chronically altered mechanical forces; if for some pathological 
reason, the remodeling of the fibers introduces some disparity in the various 
directions in the stiffness of the fibers, then it may happen that Ei 7^ E2 
and that some "dangerous" mechanical behavior develops. However, from a 
mathematical point of view the solutions of the EVPs suggest that the artery 
model is much more stable than the standard reinforcing model, due to the 
presence of exponential terms in the determining equations. 



5 Concluding remarks 

We extended the results of Merodio et al. |7j from transverse isotropy to 
orthotropy. The most important finding is that orthotropic materials may 
develop singular solutions only if there is a significant difference between the 
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mechanical stiffnesses of the two famihes of fibers. We quantified this result 
rigorously for the standard reinforcing model, where a necessary condition 
for the formation of singular solutions is that > 2/3, which means that 
one family of fibers must be at least 9.9 times stiffer than the other family. 



When we consider the arterial strain-energy density (4.1), analytical re- 
sults are no longer possible, but the methodology used to study the standard 
reinforcing material is still applicable. In this case a huge difference between 
El and E2 is necessary to possibly introduce a singularity. However, if the 
fibers are mechanically equivalent, as is usual for biological soft tissues, then 
singular solutions are avoided altogether. Therefore biological networks, such 
as the coUageneous structure of arterial walls, are the right structure to pre- 
vent the formation of the singularities described here. 

From a theoretical point of view, our results demonstrate the complexity 
of finite anisotropic elasticity and deliver some exact solutions, which are 
scarce in the literature on finite inhomogeneous deformations of orthotropic 
materials. 

It is important to note that we have barely scratched the surface of the 
collection of problems associated with the rectilinear shear of solids reinforced 
by two families of parallel fibers. Primo, we relied on strong — and perhaps, 
reductive — constitutive assumptions, namely that the strain energy density 
can be split into the sum of an isotropic part and an anisotropic part, and 
that this latter part is also the sum of two parts, each depending on only one 
anisotropic invariant. Although there is now a good body of experimental 



data supporting the adequacy of the standard reinforcing model (3.1) and 



of the biomechanics arterial model (4.1), the importance or insignificance of 



other constitutive arguments must also be evaluated, such as the role played 
by other invariants [JJj or by the angular distribution of fiber directions 
[131 [HI [IS]- Secondo, we limited our study to a shear occurring along the 
bissectrix of the two families of parallel fibers, and did not study the influence 
of other orientations. Intuitively, it is expected that this is the direction where 
the coupled reinforcing effect of the fibers is at its strongest. Nevertheless we 
were able to show that if one family of fibers is much stiffer than the other 
for the standard reinforcing model, then singularities might develop in the 
thickness of the clamped slab, in the form of discontinuities in the shear or 
in the strain gradient. 
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